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We show that full inverse covariance weighting can be naturally incorporated into modal estimation 
methods making them optimal for the CMB power spectrum, bispectrum and trispectrum, as well 
as in other 3D applications. Modal estimation methods are highly efficient requiring inversion of 
' only an rimax x rimax covariance matrix, where Timax is the limited number of modes needed to 

describe both the theoretical model and the noise which affects it (in contrast to the full /^g^^ ^ 
, lma,x inverse covariance weighting). For an isotropic modal bispectrum estimator implemented at a 

' WMAP resolution, we demonstrate how the modal inverse covariance weighting improves optimality 

significantly. We describe optimal modal estimation in general terms applicable to polyspectra of 
'~i , any order in a formalism that can be efficiently extended beyond the isotropic case. We also discuss 

further improvements by characterising and preconditioning the covariance matrix using a single 
' bispectrum variance shape induced by the noise and mask. We use the example of anisotropic noise 

in bispectrum estimation to illustrate these methods and demonstrate the necessity of subtracting 
a linear estimator term for optimality, removing cross-terms in the variance. We briefly discuss 
applications to the CMB power spectrum and trispectrum, as well as large-scale structure. 
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INTRODUCTION 



A separable eigenmode approach to CMB bispectrum and trispectrum estimation has been developed 
and implemented [1, 2], providing a reconstruction of the WMAP CMB bispectrum and a wide variety of 
^ . constraints on theoretical bispectra [3] and trispectra [4]. In this paper, we demonstrate how this modal 
Ci I methodology can be made fully optimal by incorporating inverse covariance weighting in the presence of 
anisotropic noise and other systematic effects. It can be efficiently implemented provided the eigenmode 
CN ': expansion utilised has sufficient resolution to describe, first, the theoretical polyspectrum shape being inves- 
tigated and, secondly, the noise, mask and other systematic contributions that are correlated with it. The 
Qv^ i generality and simplicity of this modal covariance inversion, combined with the efficiency of the separable 
' modal estimator, opens up the possibility of fully optimal analysis of high resolution CMB experiments such 
■ as the Planck satellite, as well as higher order correlators in the three-dimensional distribution of large-scale 
structure. 

] Full inverse covariance weighting for CMB power spectrum analysis has a long history going back to the 
^ ■ study of the COBE data using direct Cholesky decomposition with 0(/max) operations (see, for example, 
'^l [5, 6]). This was not feasible for higher resolution experiments such as WMAP, but could be improved 
• 1— I ' to ©(/^ax) with conjugate gradient methods and iterative approximations [7]. The importance of inverse 
^ covariance weighting for the CMB bispectrum has also been discussed for some time (see, for example, 
, [8]) but was only tackled at WMAP resolution more recently using multigrid preconditioning and conjugate 
gradient methods [9, 10]. Achieving full inversion at Planck resolution remains formidable goal for a realistic 
noise model with non-trivial correlations (see, for example, ref. [11]). The significant advance of the present 
work is that we show that it is not necessary to invert the full covariance matrix to obtain an optimal 
variance. Instead, we invert a much smaller modal covariance matrix from a subspace with sufficient 
degrees of freedom to describe the theoretical model, as well as the masking and noise modes that actually 
contaminate it. For example, we will investigate an isotropic bispectrum bi^i^i^ affected by anisotropic 
instrument noise or a cut-sky mask, discussing the projection of the pixel-variance noise map, variance 
and bispectrum which can affect the estimation. We have previously demonstrated, in the isotropic case, 
that a highly efficient and general modal bispectrum estimator can be implemented projecting the CMB 
data into a much smaller subspace spanned by the relevant bispectrum modes [1, 2]. What is remarkable 
is the small number of modes required to characterise almost all popular theoretical bispectrum models 
"-max = 0(10 — 100) [12]. For this isotropic subspace, we then noted that the required mode number is not 
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noticeably increased when characterising the relevant noise and mask effects [2]. Our expectation, therefore, 
when seeking an isotropic theoretical polyspectra estimator on the full anisotropic space that is that we will 
require few additional modes to incorporate relevant noise and masking effects (this will be the subject of 
a future publication). Returning to the full isotropic implementation in a WMAP-realistic context, we are 
able to significantly improve the optimality of the variance using a much simpler inverse modal covariance 
weighting, getting about 30% closer to that predicted by a Fisher matrix analysis than previous isotropic 
bispectrum estimators. 

Finally we comment on other modal estimators such as wavelets [13], bins [14] and necdlets [15] which 
differ in there choice of initial basis but are otherwise identical. Important distinctions remain, including 
the ability to analytically calculate modes, the number of modes required for convergence and the ability to 
project between primordial and CMB(or LSS)-spectra. Early versions of the wavelet method[13], claimed 
optimality while neglecting to incorporate the linear term in the estimator which is required to minimise 
the variance [16], though this has now been incorporated [17]. We demonstrate the crucial role of this linear 
term for optimal estimation in the realistic situation where rotational invariance is broken. 

The plan of the paper is as follows. We will first outline a general polyspcctrum estimation methodology 
describing projections from the full polyspcctrum space into the modal subspace of interest. In this context 
we show how the full covariance matrix is projected down into a much smaller modal covariance matrix which 
can be easily inverted. Next we discuss the modal bispectrum estimator in detail, deriving the relationship 
between the full and modal covariance matrices. WMAP results arc presented comparing analytic optimal 
Fisher matrix forecasts with the variance obtained from Gaussian simulations for the general isotropic modal 
estimator. We calculate the bispectrum variance from anisotropic noise and we note that identifying the 
primary bispectrum noise shape contributing to this variance can be used to efficiently precondition the 
covariance matrix for inversion. Finally, we give the modal estimation components required for optimal 
CMB power spectrum and trispectrum analysis and briefly discuss large-scale structure correlators. 



The CMB polyspcctrum (op) of degree p for a given non-Gaussian model is the expectation value of the 
product of p temperature multipole aim's obtained from an ensemble of universes, that is. 



where the index p represents the 2p degrees of freedom p = {Ii,mi,l2,m2, ■■■,lp,Tnp} which cover the full 
domain V of possible polyspectra (up to a given resolution / < /max)- Here we assume (a^) is the connected 
part of the correlator we seek. Most primordial theories predict an isotropic polyspcctrum which resides 
on a smaller d-dimensional subspace Vi, after integrating out (or summing over) the irrelevant anisotropic 
modes. The isotropic has an equivalent counterpart in the full space V found by multiplying by Qii± 



where p = {£,£-^} with the £ index covering the d isotropic degrees of freedom and the remaining 
2p—d anisotropic degrees. For concreteness, we note that the power spectrum Ci has d = 1 with Q££± = 
J dQxyiii,Lii^)^i,m;(^) ~ ^hh^mimii wc work with aimfi'im^ whilc for the bispectrum we have d = S with 
the Gaunt integral gei± = Qlll^ln^rn-, = ! d^xyhni,{^)Yi^m2{^)yhmA^)^ ^ ^ {hMM}, ^"^ ^ {m.i,m2,m3} 
(and similar for higher 'nondiagonal' polyspectra). 

When investigating CMB polyspectra we are typically seeking to match a statistically isotropic prediction 
Op. to the CMB data, with an appropriate weighting to reflect our knowledge of the signal-to-noise for the 
specific experiment. The estimator S determines the goodness of fit between these through the inner product 



II. MODAL CMB POLYSPECTRA ESTIMATION 



(1) 



(ap) = Gee-^ {at) 



(2) 




(3) 
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where the £„L ^ C*; JL ]>^>---C, 1, v , represents the inverse covariance across all p multipoles. Here, 
ttp" represents the 'linear term' used to subtract out known systematic effects; for the power spectrum, 
jjiin _ Ci^mi,l2,m2 ^^'^ bispectrum a''^"- = SCi^mi,l2,m20'l3m3- In practice, the linear term is calculated 

from Monte Carlo simulations of Gaussian maps with the correct experimental effects (mask, beam, noise) 
included, see [2]. We also note that (Op) is understood to be the theoretical prediction modified to include 
the relevant experimental effects. 
If we define the weighted vectors A, B and the matrix C 



.lin 

•^P ~ Ft— ' ~ 1 ^ r-i T^T- ' ^PP ~ I ' V^/ 



then the estimator can be written in matrix notation as 

X^C^^B 

Modal decompositions have proven to be very useful for estimating general CMB bispectra and trispectra. 
Suppose we have a basis TZnp where n is the mode number. We will take it to be orthonormal so 

(6) 

p 

or in matrix notation TZTZ^ = I. We can then represent our theoretical polyspectra in this basis using mode 
coefficients a defined by 

n 

The a can then be calculated easily via a = TZA. 

Most theoretical cosmological models have rotational invariance at their foundation, so the polyspectra 
A^^ they predict are generically isotropic. For this reason, our starting point must be to include sufficient 
isotropic mode functions in our basis to accurately characterise them. Most theoretical bispectra A^^ have 
smooth or simple forms and so for a suitable choice of TZ they can in fact be represented by surprisingly 
few of these isotropic modes (for example, most scale-invariant bispectrum models require only 15 modes 
at WMAP resolution [12]). This motivates us to work with a truncated basis containing rimax niodes which 
span the space of interest which is determined by the class of theoretical models under consideration modified 
by experimental effects, i.e. masking and noise. For an isotropic bispectrum estimator, we have shown that 
the experimental effects of the mask and noise can be incorporated without increasing Umax substantially 
[2]. For anisotropic contributions from experimental effects which yield a partially off-diagonal A, we need 
to supplement our basis functions with additional anisotropic modes. In line with the remarkable data 
compression possible in the isotropic case, we anticipate only increasing Umax by multiplying by a factor of 
a few provided the extra anisotropic TZ are carefully chosen. 

Given that rimax is so small relative to the vast number 2p of multipole modes, the TZ become highly 
rectangular rimaxX^P matrices. We can now define a projection operator from the full space V to the 
subspace V-p 

V = TZ^TZ (8) 

Here, we assume we have truncated our basis such that our theory and distorting experimental effects 
remain accurately described, ie. T'A = A. Using this projection we define the modal counterparts of the 
data on the same subspace V-p, 

P = TZB — >VB = TZ^P , (9) 
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as well as the covariance 



The projected covariance is related to the covariance of /3 by p! as we now show. 

{l30^)=TZ{BB'^)n^ 



^Ci^...Ci^Ci'^...Ci^^ 



(10) 

(11) 

(12) 



As we know the aim are very close to Gaussian we can expand (apttp/) into the product {aimO'l'm') A^lmO-lm) 
and {ai'm'O'i'm') ■ In general there will be p\ ways to have all pairs consist of one primed and one unprimed 
aim plus {2p)\/2p\ — p\ combinations with at least one pairing between two unprimed and two primed a;^. 
The linear term is always chosen such that it removes all contributions of the second type. As a result we 
have (BB'^) = p\C and so 



c 



p\ 



Using these quantities we can define a new estimator 



S = 



{nAfnc-^Ti^nB 

AVC-'^VB 



A^VC-^VA 



(13) 



(14) 



(15) 



where we have used {TZCVT')-^ = TZC-^TZ^ + Z_l where is an arbitrary term orthogonal to the projec- 
tion which we can ignore. Strictly C is non-invertible, formally because it contains duplicate data due to 
permutations, but also because of loss of information from a cut-sky mask. However, this problem can be 
circumvented. For example, we can restrict p to count only combinations where l\ < I2 ^ h, with a suitable 
restriction on mj when the Zj are equal, and weight them appropriately to remove this complication. Once 
it is inverted we can then reinsert duplicates to return to the previous ordering for p. Given the estimator 
expression (14), it can be easily seen that its variance is given by ct^ oc {a^C,~^a)~^ since C oc {PP'^)- 

We now wish to determine how closely the standard estimator (5) and the modal estimator (14) coincide. 
In the idealised case of homogeneous noise and no mask the covariance matrix is diagonal and so C = I 
and they are clearly identical. Also in the limit of completeness rimax^ 2p, we have V ^ I and again the 
estimators are equivalent. In a realistic context with inhomogeneous noise, a mask and a truncated set of 
modes the difference can be quantified by considering 



A 



A\\ 




B 



B\\ 
B± 



(16) 



where we define the vector = VX to be the component of X in the subspace V-p and the vector 
Afj_ = (/ — 'P)X to be orthogonal to it and, similarly, for the matrix we can decompose it into TWy = VMV, 
M-± = (/ — V)M{I — V) and A^x = VM.{I — V). In this case, we can represent the standard estimator 
(5) in the form. 



A\\ [C-'B\\ +C-/B 



(17) 
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whereas our new modal estimator (14) is 

The difference between the estimators is the inverse covariance projection of the perpendicular part of the 
data onto the parallel space. We should however pause to note that the normalisation in both cases is 
identical. As the variance we expect for the estimator is the reciprocal of the normalisation (with a factor 
for permutations), we can see the effect of the cross-projection by comparing the variance we obtain from 
simulations with this prediction. If they differ we must deduce that our covariance weighting is sub-optimal 
and the cross-term is important in which case we must increase Umax and expand our basis to include the 
relevant parts. 

It should be noted that as the covariance matrix is symmetric positive definite it is possible to find its 
Cholesky decomposition, 

C = AF. (19) 

where A is lower triangular. We can then absorb these factors into our a and P, 

a' = 'X-^a P' = A- V . (20) 

This is equivalent to re-orthonormalising our basis with respect to the full covariance matrix, and hence 
the /3 are now uncorrclatcd. This is obviously a desirable property and has several applications which we 
explore elsewhere in future publications. The most obvious is that the estimator takes a very simple form 

^ = V7 (21) 

as by construction = I. 

The computational advantages of the modal estimator (14) are considerable, especially when we note 
that direct evaluation of the full optimal estimator (5) is extremely challenging at high resolution. This is 
because, first, it represents a huge sum over /max quantities and, secondly, the inverse covariance weighting 
requires a huge I'^^g^^ x /^^j^j^ matrix to be calculated and then inverted. As we noted in the introduction, 
there do exist methods which make this tractable in specific cases, however for large /max and general noise 
models at Planck resolution this has not yet been achieved. Calculation of the mode coefficients a, (5 may 
appear to be troublesome, but it is made highly efficient using a separable form which has been discussed 
at length previously [2]. Instead, here, we focus on the improvements offered by the modal approach for 
inverse covariance weighting, reducing this to the inversion of a small nmax x 'T-max matrix. 



III. CMB BISPECTRUM ESTIMATION 
A. Inverse covEiriance weighting with the modal bispectrum estimator 

Using the CMB bispectrum, we will now demonstrate a concrete example and implementation of modal 
polyspcctra estimation with inverse covariance weighting. The estimator for the bispectrum takes the 
general form 

hmi 

where B^^jf^^^^ = {ai 

irni(^l2m2^l3m3) is the full bispectrum, modified from being purely isotropic by exper- 
imental effects, and N is the appropriate normalisation. This estimator has been shown to be optimal in 
the ideal case without mask or noise in the absence of a linear term [18]. The linear term was proposed in 



6 



[16] to minimise the variance in a realistic setting when rotational invariance is broken. To construct the 
modal bispectrum estimator, the quantities A, B, C in (4) take the form 

D h h 

A = ^ItiIIl2^ (23) 

Q 0'hm\^l2m2^hm3 ~ 3 Czirrai,f2m2^i3m3 f24) 

Q Clj^mi,l[m[Cl^rn2,l'2m'2Cl^rn3,l'^m'3 (25) 

^Ci^Ci.^ Ci.^ Ci'^ Ci'^ Ci'^ 

For simplicity and to relate the discussion more easily to previous work, let us restrict attention to an 
isotropic subspace, though the expressions derived below can be converted for the general anisotropic case. 
We shall take the usual approximation that the primary contribution to the isotropic bispectrum is the 
predicted theoretical model itself, that is, 

A = C '2 '3 ^hkh (Ofi\ 

ymimama I ^ ' 

where bij^i,^ is the predicted reduced bispectrum. Here we note that we are always working with quantities 
that include the relevant experimental effects. For example, as we arc only considering the isotropic case 
we can approximate bij^i^ ^ f skyKk2hJ}ttX ^"^"^ ~ ^^^y {bfCf"'"' + Ni) where k is the beam window 
function, Ni is the noise power spectrum and fsky is the fraction of the sky remaining after masking. 
An isotropic set of basis functions can be represented in the full anisotropic bispectrum space by 

-jphkh — ^rnim2m3 p /r,r-\ 
''^mim2m3 ~ ^hhh ^ \^') 

li I2 I3 

where vi is a weight function chosen to mimic the scaling of Q, the Gaunt integral. The simple orthogonality 
condition TZTZ'^ = / in the full space V, yields a weighted orthogonality condition on R in the projected 
isotropic subspace V-p, 

r _ \^ \Vm1m2m3) T> r> _ \^ lihh T> R (no\ 

f-^ vf vf vf ^ vf vf vf 

hrrii '1 *2 (3 ti '2 13 

where hi^i^i^ is a geometric factor defined by 



(2/i + l)(2/2 + l)(2/3 + l) (h h h 



Ku3 = V ^' ^ ) (2^) 

and decomposing the A,B,C into mode space yields the a„, Pn coefficients and the covariance matrix (^nn' 
respectively, that is. 



Olr,. = 



^ (2Zi + 1)(2;2 + l){2h + 1) f h h h V VhVl2VhKl2l3 p (or,. 



o ^miniama '*(imi<JZ2"T-2^'3m-3 Cl^rnx,l2m2^l3m3 /oi\ 

"ri= y, — /„ ^ ^ J^hhh (-^i) 



hirii 



VhVl.Vi^ ^Cl.Cifii 



nhhh (j'1'2'3 c 

i \^ ^ mim2m3^ jn'im'2m'3 p '~^limi,l['rii['~- Ijiii^-I'ji'i'j I i''>'3,l'3m'3 fo'y\ 
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Calculation of a„, /3„ is straightforward using methods described in detail in ref. [2] and, as we have chosen 
the linear term correctly to remove cross terms in the variance, = 1/6 (/3/3^) as we demonstrate below 

Umd[m[ '1 '1 '2 '3 

imi Oil2m2 [m{ 0'l'2m'2 sms^^l'^my ^^Iimi,l2m2 
m3 '^I'sm's 

) + 9Ci 

imi,l2m2 Ci' m', JLmU (a/3 ) + •••] Rnlil'2l'^ 

nhhh C C C 

„ V ^ "m.im.2'Ti3=^m', rriomo „ ^Iirnil\rn',^l2'm2,l'',m'„^l3m3,l'^m'', „ , n • 1 1 

= 6 > -Rnhhh ^ ^ -Rn'l{l'2l!, + higher Order 

= 67^C7^^ = 6C„„' . (35) 
Our inverse covariance weighted estimator (18) becomes 



(36) 



This result applies equally to the general anisotropic case where we have supplemented with basis functions 
Tlmlmlm3 which Cannot be simply represented in terms of the Gaunt integral; the derivation above is 
unaffected. Using (36) , we note that the variance cr^ of the optimal modal estimator (36) becomes simply 



2 Cnn' i^n' Pp' ) Cpp' 



<y = — rri 7^ = — — 1 • (37) 



B. Bivector representation of covariance matrix 



Further efficiencies are made possible by identifying the underlying bispectrum 'noise' and 'mask' shapes 
that contribute to the covariance matrix. We can identify such a shape because the noise is correlated 
in harmonic space. Empirically in our previous isotropic WMAP analysis we have found that the noise 
and mask contribute in tandem to produce a single shape vector P = {^n} which we could identify [3]. 
This noise/mask vector is significantly correlated with local non-Gaussianity which explains the difficulty 
in obtaining optimal constraints on this model. The dominant effect of /3 leads to the remarkable result 
that we can accurately represent the isotropic covariance matrix simply as the identity matrix perturbed 
by this bivector, that is, 

Cnn' ~ ^nn' + PnPn' ■ (38) 

For this, there is a simple inverse covariance 

- ^nn' - , Where |^|^ ^ . (39) 

1 + \pr ~1 
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For the isotropic estimator (36) in the presence of this noise/mask bivector the variance from (37) becomes 



\a\ 



1 



1 + 



1 + 



I^Pcos^e \ 



|ar V 1 + |;3|2 (1 - cos2 0) I ' 
where we have described the degree of correlation between the two vector directions through the cosine 



(40) 



a 



•/3= ^an^n = cos 6*. 

1=1 



(41) 



Prom (40), we see that the variance is made up of two components. First, there is the ideal or 'optimal' 
variance that can be achieved in the presence of the underlying Gaussian CMB signal and homogenous 

instrument noise, given the degree of sky coverage and beam resolution (incorporated into the signal-to-noise 
definition of a). Secondly, there is the contribution from the inhomogeneous noise and the anisotropic mask 
which arises if the noise/mask vector /3 is closely aligned with the model a. under investigation. Significantly, 
this second contribution vanishes if the model being studied has little correlation with j3, in which case near- 
optimal variance can be approached even in this isotropic case. In principle, in the full anisotropic case, 
the general inverse covariance weighting will approach optimality given sufficient modes for any model. 

Another key point in a practical implementation is that the single noise/mask vector /3 can be estimated 
accurately from relatively few simulations. Using (3 with the covariance matrix in the form (38) allows 
accurate and efficient inversion because we only need to determine p. This can be achieved with ngim ^ 1000 
simulations, while direct averaging of the covariance matrix takes ngim ~ 10000 to converge to percent levels 
(making it reliably invertible). We remark that early studies of these contributions in a WMAP-realistic 
context showed that the local shape is strongly correlated with the spurious contributions from the mask 
and inhomogenous noise, with the noise/mask vector /3 showing an 70% correlation with the local a^°^ 
direction. Separately, there was about an 70% local-mask correlation and a 47% local-noise correlation. 
In a forthcoming publication, we shall discuss the more general case where several 'noise' vectors can be 
obtained using a principal component analysis (PC A), separating the noise, mask and other contaminant 
vectors from model directions of theoretical interest. 

Finally, we note that for the simple bivector form of the inverse covariance matrix (39), the Cholesky 
decomposition of C can also be expressed analytically. Here, we have C, = with A the lower triangular 
matrix: 



1 + 

















\ 



(42) 



Rotation with this matrix A using (20), offers an orthonormal frame in which the bispectrum estimator 
(36) becomes simply S = a' • that is, the effects of the noise and mask are incorporated into the 

definitions of a' and fi' from the outset. 



C. Isotropic inverse covariance implementation 

The practical implementation of the inverse covariance weighting in (18) has been incorporated in the 
isotropic case for the modal bispectrum estimator. Bispectrum constraints for a wide variety of models 
using WMAP data have been presented elsewhere [3]. Here, our purpose is to focus on the improvement 
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in optimality for the variance A/nl and we only present results for local, equilateral and orthogonal non- 
Gaussian models. We compare the variance achieved for the modal isotropic estimator (which agrees well 
with previous isotropic KSW estimators) to show that it can be improved with modal inverse covariance 
weighting in this isotropic case. Of course, we do not achieve full optimality, which can be seen by comparing 
to the Fisher matrix forecast, because relevant anisotropic effects deriving chiefly from the mask cannot be 
described using only an isotropic basis. 

For the purposes of direct comparison, we present the variance AFnl using the integrated bispectrum 
Fnl defined in ref. [3], that is, the total bispectrum signal-to-noisc normalised relative to the local model 
(for which = -Fnl)- For an arbitrary predicted theoretical bispectrum bf^^^i^^, we use the estimator (22) 
to define 

F'^l = £, with N'' = N,^Nf^^)r\ (43) 
where the universal normalization N has been defined using 

with jqJ; being exactly the same quantity except defined for the /nl = 1 local model. The quantity 
Fnl provides a much more uniform variance between different models, though we also present results with 
the usual equilateral normalisation f^^^ ■ 



Estimation type 


Local AFnl 


Equilateral AFnl (A/X") 


Orthogonal AFnl (A/g-L*") 


Standard isotropic estimator 
Weighted isotropic estimator 
Fisher matrix 'optimal' 


29.5 
27.6 
22.9 


24.1 (129.0) 

23.8 (127.4) 

22.9 (122.5) 


26.0 (107.9) 
25.0 (103.8) 
22.9 (94.9) 



Table I: Results for the variance from an implementation of the isotropic modal estimator in a WMAP realistic context. These 
results were obtained from 1000 Gaussian maps at imax = 500 using the KQ75 mask and the WMAP inhomogeneous noise 
model. The mask and noise degrades the variance relative to the Fisher matrix forecast, but this is significantly improved using 
modal inverse covariance weighting, even for this isotropic estimator. 

Results for the variance AFnl from local and equilateral models are listed in Table I. This work conforms 
exactly to the WMAP-realistic methods described in ref. [3], using a KQ75 mask, WMAP beams and 
inhomogeneous noise model (see below) and with a resolution limit of Zmax = 500. The results were obtained 
using 144000 Gaussian maps from a modal bispectrum pipeline previously confirmed to be unbiased and 
which reproduces the correct /nl for simulated non-Gaussian maps [2]. Previous results for ideal full-sky 
Gaussian simulations (without inhomogenous noise or mask) demonstrate that the estimator is correctly 
normalised, confirming that optimality variance in the absence of experimental effects. For the local model 
in a WMAP-realistic context, the variance from the modal estimator is suboptimal by 30% dropping from 
the Fisher matrix forecast AFnl = 22.9 to a realised value AFnl = 29.5. This is consistent with the value 
obtained from KSW estimators previously at the same /max = 500 (e.g. AFnl ~ 26 in ref. [19]). This strong 
degradation for the local model occurs because the noise/mask vector (38)) correlates strongly with the 
local shape, increasing the variance as discussed previously (sec (40)). The addition of inverse covariance 
weighting for the local model, improves the variance to AFnl = 27.6 which takes us 30% closer to the 
projected optimal variance, demonstrating the value of incorporating this weighting even in the purely 
isotropic case. In contrast, nearly optimal results can always be obtained for the equilateral model and the 
orthogonal model because these are much less correlated with the noise/mask vector /5. For the equilateral 
model, in particular, an isotropic estimator will be capable of achieving near-optimal results. 



IV. OTHER POLYSPECTRA AND 3D APPLICATIONS 



This methodology can be extended to other correlators. We begin with the trispcctrum estimator pre- 
sented in [20], for which WMAP constraints were obtained for several models in [4]. We shall restrict 
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ourselves to the diagonal free trispectrum for simplicity, although obtaining expressions for the general case 
is straight forward. We build the optimal modal estimator our of the quantities A, B, C in (4) defined by 

nhkhU f, , , , 

A ^mim2tn3'Tt4 ^(lt2t3'4 (45) 

Q 0'limi^l2m2^l3m3^Um4, ^ Ci-^„i-^^i2m2^l3m3(^Um4, ^ Climi,l2m2 ^hmsjUmi f46) 

Ci Ci 

^yci^Ci^Ci^ Ci^ Ci'^ Ci'^ Ci'^ Ci'^ 

(^Ill2l3^4 

-r> _ ^miTn2m3m4 p /^o^ 

where vi is a suitable weight function and the isotropising factor G^^f^[\^^^^ is 

LM 

Using efficient modal methods calculation of the trispectrum estimator with inverse covariance weighting is 
no more challenging than in the bispectrum. 

Inverse covariance weighting for power spectrum estimation for high resolution experiments like Planck. 
In this case we would be constraining small deviations from a canonical model with a given Q . Considering 
an isotropic subspace, the building blocks for the model estimator are 

<^il«2 '^77117712 (Ql ~ Ql ) 

^ = 7=^ (50) 

Q _ '^hmi0-i^m2 ~ ^hmi,l2m2 /^^s 
p C'lirni,l[m[Gl2m2,l'2m'2 



(52) 



where Ci is understood to be the canonical model including the predicted deviation. Also we now consider 
Qi77ii,i277i2 — {^hmiO'l^m^ ■ '^'^^ corresponding isotropic basis functions are 

n = ^hh^I^R^i^ , (53) 

which must be supplemented to incorporate anisotropic experimental effects such as cut-sky masking. We 
can exploit the tractability of the modal inverse covariance weighting in several ways. First we could 
consider constraining specific models with specific signatures in the power spectrum with a particular an, 
such as oscillatory models. Secondly we could use the method for full power spectrum estimation using 
the recovered P to iteratively improve the best fit theoretical model. This would entail linking to Monte 
Carlo simulations which explore cosmological parameter space. Finally, we need not restrict ourselves to 
isotropic models and instead using additional anisotropic basis functions, we can to constrain anisotropy in 
the CMB. 

For large-scale structure and other 3D polyspcctra estimation, this approach has an entirely analogous 
implementation discussed in ref. [21]. We simply replace our multipoles with the density perturbation so 
now we have 

(ap) ^ (5k/k2-<^kJ (54) 
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where now the vector p represents p = (ki, k2, k„) whieh cover the full 3n-dimensional domain of 
possible polyspectra. Here momentum conservation and isotropy imply we can make the reduction to the 
subspace p = {ki,k2, ■■kn) and the two are related by 

ttp = (27r)^5(ki + k2 + ... + k„)op, (55) 

The rest of the modal method then follows analogously with two main caveats. The first being that the 
growth of structure is non-linear so even with gaussian initial conditions we produce measurable polyspectra. 
The form of these can be calculated from A^-body simulations and we can seek deviations from these due to 
primordial and other effects. Secondly there is a non-linear relationship between the late time and primordial 
polyspectra and so a Monte Carlo approach is required to find the parameter values which provide the best 
fit to the data. Error bars can then be obtained from simulations with the chosen primordial spectrum. 
However, the inverse covariance matrix can still be calculated from Gaussian simulations, so it does not 
need to be recalculated for differing levels of primordial non-Gaussianity. 

V. CONCLUSION 

We have shown that full inverse covariance weighting can be naturally incorporated in the modal polyspec- 
tra estimator to achieve optimality. The resulting modal covariance matrix can be efficiently calculated and 
inverted. This greatly reduces the computational cost and memory requirements of polyspectra estimation, 
making optimal analysis of high resolution experiments like Planck possible. We emphasise that the pri- 
mary efficiencies discussed here arise because we have effectively projected out all but the modal degrees 
of freedom needed to characterise the model (or models) under investigation, while retaining an adequate 
description of the noise and systematics with which it is correlated. 

We have only briefly discussed the many applications of this work. These include improved WMAP7 
constraints on a wide variety of theoretical models, together with a model-independent constraint on the 
overall magnitude of the bispectrum in ref . [3] . A principal component analysis can be used to identify the 
noise, mask and other shapes for preconditioning the inverse covariance weighting. A forthcoming paper will 
use eigenmodes to characterise and marginalise over foregrounds and other contaminants in preparation for 
the CMB analysis of Planck data. The analysis of higher order correlators in three dimensions represents 
a particular challenge as many billion galax:y surveys become available. The 3D modal estimator [21] 
allows for efficient polyspectra estimation and for the setting of arbitrary non-Gaussian initial conditions. 
Combined with inverse covariance weighting, this methodology may help unravel the many systematic and 
evolutionary effects present in large-scale structure observations. 
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Appendix A: Anisotropic noise and bispectrum optimality 

Anisotropic instrument noise provides an explicit example which illustrates some of the points raised 
above, notably the nature of its contribution to the variance the importance of the linear term in (5) for 
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optimal estimation. WMAP non-Gaussian analysis to date has used a simple instrument noise model which 
assumes no correlations between pixels. If n/^ is the harmonic transform of the pixel-variance map, then 
the covariance of the noise contribution to the CMB temperature is given by 



47r 



(Al) 



(See, for example, the discussion of covariance under rotations in [22].) The covariance of the cut-sky mask 
also can be represented in this form by taking the masked pixels to have a very large noise so that the 
covariance weighting neglects that section of the map. We shall now consider the contributions to the 

variance of the estimator both with and without the linear term. 

We assume CMB signal plus instrument noise can be represented as ai„i + aJ^^ with respective covariances 



hma.m. = {ahm,ai,m,) = {-l}'"' QJij^^^.^^ and C;^^^ ^^^^ = {af'^^^afl^J. Again restricting attention 
to the isotropic subspaee by expanding using isotropic modes (27), the modal covariance matrix with the 
inhomogeneous noise takes the form 



Cnn' — 6 ^nlM3Qrnim2m3^limi,l'im'iCi2m2,l'2m'fihm3 
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where Cj^ represents the isotropic power from the pixel-variance map. 



C 



E 



21+1 



'^Im'^l—m ; 



^fihh '^^ isotropic part of its bispectrum, 



mi 



(A2) 



(A3) 



(A4) 



Standard identities were used to convert summed products of four and five Wigner-3j symbols into the 
Racah-6j and Racah-9j symbols respectively. 

The expression (A2), which yields that the isotropic covariance from anisotropic noise (Al), depends 
solely on the two- and three-point correlators of the pixel-variance map (i.e. the isotropic functions Cj^ 
and 6^;^; ). The reduction to isotropy means that the product form of the full covariance matrix C can be 
reduced from a Z^gx ^ ^max matrix to x /^gx but this is still intractable. Here, however, we have a 
modal covariance matrix of only nmax x 'i-max which can be easily inverted. 

Now consider the contribution from inhomogenous noise to the covariance matrix in the absence of a 
linear term in the estimator (31) but without the linear term. In this case, the full covariance has the 
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following additional cross-terms 



"(2/i + l)(2/i + l 



47r 



+ 



2h^ /j2 



(A5) 



The first line of the resulting expression corresponds to a monopole so can be ignored. The remaining two 
lines containing three cross terms do not in general vanish, but once again they contain only contributions 
from the isotropic parts of the pixel- variance map correlators Cf^ and ^^i^i^- Other effects which break 
rotational invariance in a similar manner, such as a cut-sky mask, can be treated similarly. 

It is important to note that the existence of the non-zero cross-terms (A5) means that the minimum 
variance cannot be achieved without use of the linear term as shown in (33). 
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